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We investigate excited states of the Shastry-Sutherland model using a kind of variational 
method. Starting from various trial states which include one or two triplet dimers, we 
numerically pursue the best evaluation of the energy for each set of quantum numbers. 
We present the energy difference as a function of either the coupling ratio or the mo- 
mentum and compare them with the perturbative calculations. Our data suggest that 
the helical order phase exists between the singlet dimer phase and the magnetically or- 
dered phase. In comparison with the experimental data we can estimate the intra-dimcr 
coupling J and the inter-dimer coupling J' for SrCu2 (1303)2: J'/ J — 0.65 and J = 87K. 
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1. Introduction 

One of recent fascinating topics in two-dimensional quantum spin systems is the Shastry- 
Sutherland model (SS model), ^ the model of orthogonal dimers with the intra-dimer coupling 
J (> 0) and the inter-dimer coupling J' (> 0). This model provides a nontrivial system which 
relates an exactly solvable spin-gaped model obtained in the J' — > limit to a square lattice model 
where intra-dimer interaction vanishes (J = 0). It is known that the exact ground state of the SS 
model is the direct product of the singlet dimers when J' < J' c ~ 0.68J, while for the sufficiently 
large J' /J the system is in the magnetic order phase with the Neel ground state. 

A lot of notable studies have been done on this system since Kageyama et al. 2 ^ found that 
SrCu2(B03)2 realizes the SS model and its coupling ratio J'/ J is very close to the critical value of 
the model. Active experimental works have reported additional interesting features of this matter 
such as magnetization plateaus, 2 ~ 4 ) magnetic-field dependence of the thermal conductivity, 5, 6 ) and 
soundwave anomalies. 7 ) 

Theoretical investigation was stimulated by Miyahara and Ueda 8 - 1 who pointed out the locality 
of the excited states in the SS model. Owing to the expectation that this locality should account 
for the observed magnetization plateaus, many calculations based on the perturbative methods 9-14 ) 
have been carried out. Although their results qualitatively agree with the experimental data, 15-18 ) 
the expansion parameter, estimated to be 0.6 — 0.7 for SrCu2(B03)2, seems too large to ensure the 
validity of the perturbation. Therefore non-perturbative approaches are desired in order to confirm 
discussions on the excited states. 
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In this paper we show numerical results obtained by a variational method, which we have 
recently developed and named the operator variational (OV) method, 19 ' 20 ) together with the re- 
structuring technique. 21 ' 22 ^ Our purpose is to study the lowest excited state of the SS model on 
a square lattice sketched in Fig. 1. Using the OV method we can perform systematic calculations 
without damaging the quantum symmetries of the system. The trial states, on the other hand, 
should be given by inspection. We examine some of the one-triplet states and the two-triplet states 
with the total spin S = 0, 1 and 2, paying special attentions to the symmetries of the system's 
Hamiltonian on lattices equally sized in the x and the y directions. We see that the value of the 
lowest energy, which provides an upper bound of the true value, strongly depends on the choice 
of the trial state. We therefore employ various trial states and adopt the one whose energy is the 
lowest among them. 

Our results indicate several interesting features on the dispersion relations and on the phase 
transition of the model. They are in good accordance with the perturbative results if the inter- 
dimer coupling is weak, but discrepancies are observed when we move on toward the transition 
point. It is likely that the helical order phase 23 ' 24 ) exists. Trying a comparison with experimental 
data on SrCu2(B03)2 we obtain several encouraging results. 

In section 2 we introduce the SS model and discuss its symmetry used to distinguish the states. 
Section 3 is to give a brief description on the OV method and to comment on its application to the 
SS model. Numerical results on 8 2 and 12 2 lattices are presented in section 4 and the final section 
is devoted to summary and discussions. In appendices A-E we show the trial states employed in 
our numerical study. 



2. Model and Symmetry 

In the Shastry-Sutherland model a spin is located at every site of the lattice schematically 
shown in Fig. 1, whose coordinates are denoted by either 

{2n x a ± ^, 2n y a ± ^) (A sublattice) or ((2n x + l)a ± ^, (2n y + l)a =F ^) (B sublattice) 

with n x ,n y = 0, 1, • • • ,N e ff — 1, where N e ff is the number of dimers, 2a is the unit distance 
between dimers and d is the distance between two spins of a dimer. The total number of spins N 
therefore equals to 4iV 2 jj. 

The Hamiltonian of this model is given by 



H = -J ^2 {h a {n x ,n y ) + h b (n x ,n y )} 



4 

n x ,n v =0 



+ jJ' ^2 {h 1 (n x ,n y ) + h 2 (n x ,n y ) + h 3 (n x ,n y ) + h4(n x ,n y )}, (1) 



N eff -1 
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where, cr(x,y) being the Pauli matrix at the location (x,y), 

h a (n Xl n y ) = a(2n x a + ^,2n y a + |) • a(2n x a - ^,2n y a - |), (2) 

h b (n x ,n y ) = cr((2n x + l)a + ^,(2n y + l)a-^)-<T((2n x + l)a-^,(2n y + l)a + ^), (3) 

d d 
hi{n x ,n y ) = <r(2n x a+ -,2n y a+ -) ■ 

{<r({2n x + l)a + |, [2n y + l)a - |) + cr((2n ;c + l)a - |, (2 % + l)a + |)}, (4) 

d d 
h2{n x ,n y ) = a(2n x a - -,2n y a - -) ■ 

{a((2n x - l)a + |, (2n y - l)a - ^) + <r((2n :c - l)a - ^, (2 % - l)a + ^)}, (5) 
h 3 (n x ,n y ) = a((2n x + l)a + ^, (2n y + l)a - ^) • 

{cr((2n ;E + 2)a + |, 2 % a + |) + <r((2n x + 2)a - ^ 2nj,a - ^)}, (6) 
h^(n x ,n y ) = cr((2n x + l)a - ^, (2n y + l)a + ^) • 

{(7(2^0 + ^, (2 % + 2)a + ^) + o-(2n a: a - |, (2 % + 2)a - ^)}. (7) 

Periodic boundary conditions are assumed in both x and y directions. 

Symmetric operates to commute with the Hamiltonian are T x (translation in the x direction), 
T y (translation in the y direction), U (A-B translation), V x (A-B translation in the x direction), 
V y (A-B translation in the y direction) and I± (x-y reflections), which translate cr(x,y) as follows. 

T x a{x, y)Tl = <t(x - 2a, y), T y a(x, y)T% = a(x, y - 2a), 

Ua(x,y)U^ = cr(—y — a,x + a), 
V x cr(x, y)V% = a(-x -a,y + a), V y a(x, y)V y ^ = a(x + a,-y- a), 
I+<r(x, y)!+ = (r(y, x), I-<t(x, y)I ] _ = tr(-y, -x). 

Note that 

U 4 = l, 4 = 1, I 2 _ = l, 
where 1 denotes the unit operator, and 

U = I+Vy = I-V x . 

3. OV method 



In this section we make a brief description of the OV method we proposed in ref.19, where we 
showed that this method enables us to calculate the energy for large quantum systems with less 
computer memory resources compared with the Lanczos method. 
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In the OV method, assuming that the Hamiltonian H is the sum of N g partial Hamiltonians 

9i, 

i=i 

we systematically generate candidate operators 

N a 

d(k\ = 0, fc 2 , • • • , &l) = S ^9i+k 1 9i+k 2 • --9i+k L (L = l,2,--- , L max ), 

i=i 

with some integer L max and all possible values of integers &2, ■ ■ ■ , Then we select a set 
of operators {Oj} = {Oo = I(Identity),Oi, ■ ■ ■ , O n -i} from those candidates so that Oo I ^)) 
Oi | VP), ■ ■ ■ , O n -\ | VP) includes all linearly independent states for the given L max and | VP). Thus 
we obtain an approximate basis to calculate the Hamiltonian matrix elements. Once the basis is 
determined, the orthogonalization and the diagonalization are easily carried out by conventional 
methods. 

Let us then concentrate our attention to trial states, the highly model-dependent component of 
the method. We study the dispersion in two cases. One is the case p y = p x , where we employ those 
states which have even /odd parity in the x-y reflection I + . They are one-triplet states with the 
total spin S = 1, two-nearest-neighboring-triplet states with S = 0, 1, 2 (abbreviated as nn-triplet 
states hereafter) and two-next-nearest-neighboring-triplet states (nnn-triplet states hereafter) with 

5 = 0, 1, 2. In the other case p y = 0, the trial states are eigenstates of V y , the operator to represent 
A-B translation in the y direction. Details of these trial states are given in appendices. 

4. Results 

Now we present our results obtained on 8 x 8 and 12 x 12 lattices (p x or p y being therefore 0, 
|7r, i"7r, |7r and 7r) with L max = 3. We did not find any size effect on the results with the same 
(Px, Py), where p x (p y ) = or n, calculated for these two lattice sizes. 25 ' 26 ) 

In Fig. 2 we plot, as a function of the coupling ratio J'/ J, the difference between the energy 
of the singlet-dimer state Eq and the lowest energy E for the momentum p=(0, 0), the total spin 
S = j (j = 0, 1) and the even or odd I + parity. We determine E with S = 1 as the minimum 
value among the energy calculated from the one-triplet trial states, the nn-triplet trial states and 
the nnn-triplet trial states which are explicitly described in the Appendices D and E. Let us 
introduce notations E™^, E™ n and to represent the minimum energy obtained from the 

one-, the nn- and the nnn-triplet trial states, respectively. Then the above is stated by E = 
min{E™? n , E™ in , E™£} . As for the S = states we compare the nn- and the nnn-triplet trial 
states, namely E = min{E?™ n , E™%} . 

When p=(0, 0) it turned out for S = 1 that E = E™ in and E is doubly degenerate, being 
independent on the 1+ parity. Our results on the energy difference E — Eq, which correspond to the 
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spin gap for this value of S, nicely agree with the perturbative expansion up to the fifth order of 
J'/J (the dotted line in the figure) presented by Fukumoto. 13 ) Especially for small values of J' /J, 
say for J'/J < 0.5, the agreement is excellent. When S = 0, E is obtained from the nra-triplet trial 
states for the odd I + parity in the whole range of the examined coupling ratio, while for the even 
I + parity E = E™-™ if J'/J > 0.55. Both I + parity results are also compatible with the ones based 
on the perturbation theory, 12 ). 14 ) 

A remarkable feature in Fig. 2 is a crossover between E with S = 0, odd I + parity and E with 
5 = 1. We see that the former, which is much larger than the latter when J'/J is 0.45, decreases 
more rapidly to reach Eq faster than the latter when the coupling ratio grows. This observation 
suggests the existence 27 ) of the helical order phase 23 ' 24 ) between r\ < J'/J < where r\ (r 2 ) is 
the value of the coupling ratio where E with S = 0, odd I + parity (with S = 1, even and odd I + 
parity) becomes equal to Eq. There are two reasons to believe that this crossover is not an artifact 
due to the approximations in the OV method but a real property of the SS model. One is that the 
exact calculations on small lattices (16 and 24 sites) also indicate this crossover. Another reason is 
that the S = 1 results would be less affected by the higher order corrections than the S = ones 
because they are known to be more reliable in the coupling region around the crossover. Since the 
OV method is a kind of variational method, this means that we would observe more decrease of E 
in the S = case when L max is increased. The crossover is thus expected to survive in more precise 
estimations. Quantitatively, the crossover point and the values of r\ and r 2 , which are ~ 0.65, 
~ 0.75 and ~ 0.8 in the present study, might change for the smaller in further numerical work. 

Next we study how E depends on the momentum (p,p) or (p, 0). In Figs. 3-5 we show 
dispersions of the energy difference calculated for p=(p,p), S = j (j = 0, 1, 2) and p=(p, 0), S = 1 
with the coupling ratios J'/J = 0.5 (Fig. 3), 0.6 (Fig. 4) and 0.65 (Fig. 5). Here we examine 
the one- and the nn-triplet trial states for S = 1 and p=(p, 0) , while for p=(p,p) we calculate 
E = min{E™ n , E™} {E = min{E™f n , E™ n , E™}) when S = 0, 2 (S = 1). 

Let us first direct our attention to the dispersions from the S = 1 trial states. When p=(p,p) 
we observe, at each coupling ratio stated above, that the lowest energy with even I + parity is 
^rnin ^ or an va l ues of p. With odd I + parity, on the other hand, we observe that E = E^f n for 
< p < it. As was stated above, the lowest energy for p = is irrelevant to the 1+ parity so that 
the system has doubly degenerate E(= E™ in ). This degeneracy is also found at p = it. In another 
case p=(p, 0), we find that £™J n is slightly smaller than E°^ n when p = tt while E^ n > E^ n for 
all other values of p. The S = 1 results in the figures indicate that the even I + state is almost 
dispersionless. Its value changes only ~ 3 % even when the value of J'/J is as large as 0.65. We also 
see that the energy with the odd 1+ state depends very weakly on p. Similar tendency is observed 
for the p=(p, 0) state, too. This feature is in marked contrast to the reported perturbative work, 
where the dispersion is flat only for small values of the coupling ratio. For example, the dispersion 
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by Weihong et al. 10 ^ shows rapidly growing fluctuations near J' / J ~ 0.65. 

Turning to results for 5 = state we see strong momentum dependence of E both for even 
and odd I + parity, which is the very impressive characteristic of this state. Here E = £™^ n always 
holds if the i+ parity is odd. For the even I + parity, however, it depends on the values of J' / J 
and p whether E™ in < E™%% or not. When the coupling ratio is 0.5 it turned out that E = E™ n 
whatever p is. On the other hand, we observe E%£™ underlies EJ™ n for J' / J = 0.6, p = and for 
J' I J = 0.65, p = 0orp = 7r/3. These dispersions are in qualitative agreement with those given 
by Totsuka et aZ. 12 ) for J'/J=0.5 and by Fukumoto 13 ) for J'/J=0.55, although their values are 
smaller than ours. 

Finally we comment on the 5 = 2 case, where the relation E™£ < E™ in is always observed. 28 ) 
Both the even and the odd I + parity dispersions are weakly dependent on the momentum in 
accordance with statements in refs. 12 and 13. Our results did not confirm, however, their claim 
that the energy difference is smaller than the twice of the spin gap so that the 5 = 2 state is a 
stable bound state composed of the lowest excited state. What we read from Figs. 3-5 is that 
E — Eq for 5 = 2 is located slightly above the twice of that for 5 = 1. 

Before proceeding to the next section we would like to make a comparison between our results 
and the experimental data for SrCu2(BOs)2. The energy of the lowest excited state with 5 = 1, 
namely the spin gap, can be determined in a very reliable manner both in experiments and in 
numerical works. The spin gap for this matter is established to be 3.0 meV (24.2 cm -1 , 34.7 
K) through various experiments. 4 ' 16 ' 17 ) Another reliable quantity shared by experimental and 
theoretical investigations would be the energy of the lowest excited state with 5 = 0, which is 
reported to be 30 cm -1 (3.7 meV) in the Raman scattering experiment. According to Knetter et 
a/., 14 ) the observed energy is that of the SD=— state (the odd U parity state in our terminology) 
because the SD=+ state, which has the lower energy, is forbidden in this experiment. Among 
our trial states with 5 = and p=(0,0) shown in the Appendix D, the odd U(= I-V x ) ones 

are I *rm,0,(0,0),a')> I *wm,0,(0,0),& / ) and I *wm,0,(0,0),c') wmle I *ran,0,(0,0),&')> | *„nn,0,(0,0),o') and 

I *rmn,o,(0,o),d') ^ ave the even U parity. What we observe are that the lowest (the second lowest) 
energy is -E n „, ,(o,o),a' (^nn,o,(o,o),fc') if th e coupling ratio is less than or equal to 0.8 (0.7). Thus, 
although we plot ^nn,o,(o,o),o' m Fig- 2, we should examine £? n n,o,(o,o),b' here. Requesting that the 
ratio (-E n n,o,(o,o),6' — -^o)/(-Snn,i,(o,o),c' — Eq) should be equal to the experimental value 3.7/3.0 = 
1.23, where -Enn,i,(o,o),c' is the lowest energy with 5 = 1, we fix J'/ J for SrCu2(B03)2 is 0.65 and 
J = 87K. These values are very close to those estimated by in ref.29, which are 0.635 and 85K, 
respectively. 

Table 1 lists up our results for 5 < 1 at J'/ J = 0.65 which are the lowest energy with different 
symmetries. From Table I we can further estimate, using the value J = 87K, that the energy 
difference of the odd U, even I + state (of the odd I+I- state) is about 6.8 meV (7.0 meV). These 
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might correspond to the 5 = state observed in the infrared experiment, 16 ) whose energy is 52 
cm -1 (6.4 meV), or the one in the Raman scattering, 15 ) whose energy is 56 cm -1 (6.9 meV). 
We did not find, however, any 5 = 1 state to explain the second excited state observed in the 
experiments 4 ' 17 ^ with the energy 54. 7K (4.71 meV). The reason would be that the OV method is 
a kind of variational approach where the higher excitations for each set of quantum numbers are 
difficult to calculate. 

5. Summary and discussions 

In this paper we made a report of the excited states of the SS model. We carried out numerical 
investigations using the OV method, a kind of variational method we have developed. 19 ' 20 ) Our 
results show that the dispersion of the lowest excited state with the total spin 5 = 1 is quite small, 
which indicates the locality of the state, even when the coupling ratio J' /J is as large as 0.65. It 
provides a contrast to the perturbative results, where the small dispersion is observed only for small 
J'/ J. 10 ) Another noticeable feature in our observations is that, when we start from J' /J = 0.45 
and increase this ratio keeping the momentum p=0, the energy of the 5 = state decreases faster 
than the energy of the 5 = 1 state to form a crossover near J' / J = 0.65. This crossover supports 
an existence of the helical order phase suggested in refs. 23 and 24. The upper bound of J' /J 
for the phase transition point between the orthogonal dimer phase and the helical order phase is 
0.75 in this study. We also successfully estimate the values of J' / J and J from our results and 
the experimental data, J' / J = 0.65 and J = 87K for SrCu2(B03)2, which are compatible with the 
preceding studies. 14, 29 ) 

Since it is so far difficult to perform the next order {L max = 4) calculations, which would be 
the best way to figure out the errors in our results, let us here try to discuss the reliability of 
our data (labeled by OV hereafter) by comparing some of them with the perturbative results 12-14 ) 
(PB12,PB13,PB14), with the values from the perturbative expansion on the spin gap up to the fifth 
order of J' /J 13 ' (PBF) and with the results from the exact diagonalization of 24 or 16 site cluster 
(ED24, ED16). For the lowest 5 = 1 state, we have 0.679 (OV), 0.685 (PBF), 0.6779 (ED24) and 
0.6797 (ED16) at J' /J = 0.5 while they are 0.402 (OV), 0.436 (PBF), 0.3676 (ED24) and 0.3836 
(ED16) when J' /J = 0.65. In the similar analysis on the lowest excited state with 5 = 0, we obtain 
0.93 (OV), 0.80 (PB12), 0.90 (PB13, PB14), 0.858 (ED24) and 0.817 (ED16) at J' /J = 0.5 while 
they are 0.59 (OV), 0.55 (PB13), 0.50 (PB14), 0.432 (ED24) and 0.386 (ED16) when J' /J = 0.6. 
These results suggest that our data are quite reliable when J'/ J ~ 0.5 and for larger values of the 
coupling ratio they give us reasonable upper bounds of the energy difference, especially for the spin 
gap. 

We did not observe any instability on the 5 = 1 two triplet states for large J' / J in contrast to 
the report in refs. 13 and 14. Because this is an important issue in connection with the existence 
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of the helical order phase, further study would be necessary before we draw a definite conclusion. 
Another physically intriguing phenomenon is the magnetization plateau of SrCu2 (603)2- It is so 
far an open question to explain the experimentally observed 1/8 plateau from a theoretical point 
of view. 11 ' 30 ) We hope our method is helpful to solve this problem. 

Appendix A: Ground state 

It is known that the exact ground state | ^0) f° r J' /J < 0.68 is the product of all singlet 
dimers, 

I ^0) = J I I singlet(2n x ,2n y ))- \ singlet(2n x + l,2n y + 1)), 

n x ,n y =0 

with 

I singlet(2n x ,2n y )) = 

1 ,.. d d d d 

-j={\] (2n x a + -, 2n y a + -) [ {2n x a - -, 2n y a - -)) 

d d d d 

- || [2n x a+ -,2n y a+ -) | (2n x a - -,2n y a- -))}, 

I singlet(2n x + 1, 2n y + 1)) = 
-^={|T {{2n x + l)a + ±, (2n y + l)a - ^) [ ((2n x + l)a - ±, (2n y + l)a + ^)) 

- |I (pn* + l)a + ^, (2n y + l)a - ^) ] ({2n x + l)a - ^ (2n y + l)a + ^))}. 

Appendix B: Operators which create the triplet dimer states 

We introduce the operators Tj (j = 0, ±1) which operate the dimer singlet state and create 
the triplet one, 

Tj(2n x ,2n y ) \ signlet(2n x ,2n y )) =| S z = j(2n x ,2n y )), 

fj(2n x + 1, 2n y + 1) | signlet(2n x + 1, 2n y + 1)) =| S z = j(2n x + 1, 2n y + 1)), 



where 



fi{2n x , 2n y ) = -^={-a + (2n x a + ^, 2n y a + |) + <r + {2n x a - |, 2n y a - ^)}, 
f x {2n x + 1,2^ + 1) = 

-L{-a + ((2n x + l)a + |, (2n y + l)a - |) + <7+((2n x + l)a - |, (2 % + l)a + |)}, 
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- 1 d d _ d d 

T-i(2n x ,2n y ) = -—{-a {2n x a + -, 2%a + -) + a [2n x a - -, 2n y a - -)}, 

f_ 1 {2n x + 1,2% + 1) = 

-^{-^((2n s + l)a + ^ (2% + l)o - ^) + a-((2n x + l)a - |, (2% + l)a + ^)}, 

T (2n x ,2n y ) = --{-a (2n x a + -,2n y a + -) + a (2n x a - -, 2n y a - -)}, 
Tb(2n x + 1,2% + 1) = 

-^{-**((2n* + l)a + ^, (2n y + l)a - |) + a z ((2n x + l)a - |, (2% + l)a + |)}. 

Appendix C: States in the momentum space 

We define the states in the momentum space as follows. 

I (Px,P y )A) = ^2 cxp(ip x n x + ip y n y )fi(2n x ,2n y ) | * ), 

Tlx ifly 

I (Px,P y )B) = ^2 exp(ip x n x + ip y n y )f 1 (2n x + l,2n y + 1) | * >, 



Tlx ifly 



I (p x ,p y )AR2) = ^2 exp(ip x n x + ip y n y )f 1 {2n x ,2n y )f 1 (2n x + l,2n y + 1) | * ), 

Tlx 

I {Px,P y )AL2) = ^2 exp(ip x n x + ip y n y )f 1 (2n x ,2n y )f 1 (2n x - l,2n y - 1) | * ), 

fix ifly 

| {p x ,p y )BR2) = exp(ip x n x + ip y n y )f 1 (2n x + l,2n y + l)f 1 (2n x + 2,2n y ) | 

Tlx r^y 

| (p x ,p y )BL2) = ^2 exp(ip x n x + ip y n y )f 1 (2n x + 1,2% + l)f 1 (2n x ,2n y + 2) | * ), 

Tlx ifly 

I {Px,P y )ARl} = ^2 exp(zp x n x + ip y n y ) ■ 

Tlx }Tiy 

{fi(2n a! , 2%)T (2n x + 1, 2% + 1) - T (2n x , 2 % )f 1 (2n x + 1, 2n y + 1)} | * ), 
| (p x ,p y )ALl) = ^2 exp(ip x n x + ip y n y ) ■ 

fix ifly 

{Ti(2n a .,2n w )T (2ra x - 1,2% - 1) - T (2n x , 2n y )f 1 (2n x - 1,2% - 1)} | * ), 
I {p x ,p y )BRl} = ^2 exp(ip 

x^x H~ ^Py^y) 

fix ifly 

{Ti(2n x + 1,2% + l)f (2n x + 2,2%) - f (2n x + 1,2% + l)Ti(2n x + 2,2%)} | *„}, 
I {Px,P y )BLl} = ^2 exp(ip x n x + ipyUy) • 

fix ifly 

{Ti(2n x + 1, 2% + l)f (2n x , 2% + 2) - f (2n x + 1, 2% + l)Ti(2n x , 2% + 2)} | ^ ), 
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{Ti(2n x , 2n y )f- 1 (2n x + 1, 2n y + 1) + T_i(2n x , 2n y )Ti(2n x + 1, 2n y + 1) 
-f (2n x , 2n y )f (2n x + 1, 2n y + 1)} | * ), 
I (Px,Py)AL0) = ^2 exp(ip x 

Tlx H~ ^PyTly ) ' 

{fi(2n a .,2n w )f'_i(2n a . - l,2n y - 1) + T_i(2n x , 2n y )Ti(2n x - 1,2 % - 1) 
-f (2n a: ,2nj / )fo(2n :r - l,2n y - 1)} | * ), 
I (p x ,Py)BR0) = expfe 

IT'X r^y 

{Ti(2n x + 1, 2n y + l)T_i(2n x + 2, 2^) + T_i(2n x + 1, 2n y + l)Ti(2n x + 2, 2n y ) 
-f (2n x + 1, 2 % + !)%{2n x + 2, 2n y )} | tf >, 
I (Px,P y )BL0) = ^ exp(zp x 

IT'X y^y 

{Ti(2n x + 1, 2n y + l)T_i(2n x , 2^ + 2) + T_i(2n x + 1, 2^ + \)f x (2n x ,2n y + 2) 
-T (2n x + 1, 2n y + l)T (2n x , 2 % + 2)} | * ), 

I {p x ,p y )Ax2} = ^ exp(ip x n x + ip y n y )f 1 (2n x , 2n y )f 1 (2n x + 2, 2n y ) | * ), 
I (Px,Py)Ay2) = ^ exp(ip :r n :E + ip y n y )f 1 (2n x , 2n y )f 1 (2n x ,2n y + 2) | * ), 

n x ,n y 

| (p x ,p y )Bx2) = ^ exp(ip a; n :r + ip y n y )fi(2n x + l,2n y + 1)71 (2n x + 3,2n y + 1) | VJ> ), 

f^x^y 

| (p x ,p y )By2) = ^ exp(ip x n :r + ip y n y )fi(2n x + l,2n y + 1)71(2723; + l,2n y + 3) | * ), 
I (PccPj/)^!) = ^ exp(ip x n x + ip y n y ) • 

{f 1 (2n x ,2nj / )fo(2n a: + 2, 2n y ) - T (2n x , 2n y )Ti(2n x + 2, 2n y )} | *„>, 
I (Pa;,Pj/)^yl) = X] e *v{iPxn x + ip y n y ) ■ 

{f 1 (2n x ,2n y )f {2n x ,2n y + 2) -f Q (2n x , 2n y )f 1 {2n x , 2n y + 2)} | * ), 
I (Pa;,Pj/)-Ba;l) = ^ exp(ip x n x + ip y n y ) ■ 

f^x y^y 

{Ti(2n x + 1, 2n y + l)f (2n x + 3, 2n y + 1) - f (2n x + 1, 2n y + l)Ti(2n x + 3, 2n y + 1)} | *„>, 
I {Px,Py)Byl} = ^ exp(^ x n x + ip y n y ) • 

{Ti(2n x + 1, 2n y + l)f (2n x + 1, 2n y + 3) - f (2n x + 1, 2n y + l)Ti(2n x + 1, 2n y + 3)} | * )- 
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Appendix D: Trial states in the p y = p x case 

Using the states defined in the previous section we construct trial states. We employ S z = S = 
j (j = 0, 1,2) states described below as the trial states. It should be noted that some states are 
degenerate. In order to point out the degeneracy we use the notation E m to represent the lowest 
energy obtained with the trial state | * m ). 

D.l one-triplet states (S = 1) 

The trial states with one triplet are as follows. 

I *one,i,(pj»),a) = I (p>p) A ) (I+even), 
I ^one,l,(p,p),b > = I (P,P)B) (I+odd). 

Note that £? one ,l,(p,p),a = E<me,l,(p,p),b- 

The following should also be kept in mind in order to make a comparison between the numerical 
results and the experimental data. 

I- I *orae,l,(0,0),a) = — I ^ one, 1,(0,0), a) Vc | ^ one, 1,(0,0), a) = — I *orae,l,(0,0),&)> 
I- I ^ one, 1,(0,0), b) = + I ^ one,l,(0,0),b) i Vc | * on e,l,(0,0),b) = ~~ I *one,l,(0,0),a) ■ 

D.2 two-nearest-neighboring-triplet states (S = j, j = 0,1,2) 

The trial states with two nearest neighboring triplets for p / 0, p / -it are 

I ^nn,j,(p,p),a) = I (P,P)ARj) (I+odd), 
I ^nn,j,(p,p),b) = I (.P,P)ALj) (I + odd), 
I ^ nn,j,(p,p),c 

) = | (p,p)BRj) + | (p,p)BLj) (I + odd), 
I ^nn,j,(p, P ),d) = I (p,p)BRj)- | (p,p)BLj) (I+even). 
For the momentum (0,0) they are 

I *„„j,(o,o) ) a') = I (0, 0)ARj)+ | (0, 0)ALj> + | (0, Q)BRj)+ | (0, 0)BLj) (I+odd), 
I *nnj 1 (o,o),6') = I (0, 0) ARj>+ | (0, 0)ALj) - I (0, 0)5i?j> - | (0, Q)BLj) (I+odd), 
I *„„j 1 (o,o),c') = I (0,0)Ai?j)- | (0,0)ALj) (I+odd), 
I *„„j,(o,o) ) d') = I (0,0)Si?i)- | (0,0)BLj) (I+even). 
Also for the momentum (it, 7r) they are 

I *nn 1 j,( w , 1 r),o") = I (?r, tt) Ai?j)+ | (it , ir)ALj) + i{\ (ir,ir)BRj)+ I (7T,7r)SLj)} (I+odd), 
I ^nn,j,(^),b") = I (vr,^)Ai?j)+ I (tt, ir)ALj) - i{\ (ir,n)BRj)+ \ (ir,it)BLj)} (I+odd), 

I ^nn,j,(n,n),c") = I {lT,n)ARj)- | (?T, 7r)ALj) (I+odd), 

I *nn,j,(7r,7r),d") = I (vr, ir)BRj)— \ (ir,Tr)BLj) (I+even). 
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Note that £ n nJ,(0,0),c' = #rm,.j",(0,0),d' and E nn j^^),c" = E nn,j, (%,%), d" ■ 
The following should also be kept in mind, too. 

I- I ^nn,j, (0,0), a') = _ I ^ nn,j, (0,0), a') ■> I ^ nn,j, (0,0), a') = + I *nn J, (0,0), a'} j 

^- I *nn,j,(0,0),b') = _ I *nn J,(0,0),&') > I *nn,j,(0,0),6') = _ I *rmj",(0,0),&')> 

^- I *nnj,(0,0),c') = + I ^ nn, j, (0,0), c')i Ke I *nnj,(0,0),c') = ~~ I ^ nn,j, (0,0), d') ■• 
I- I ^ nn,j, (0,0), d'} = _ I *nn,j',(0,0),d')) I *nn,j,(0,0),<2') = — I ^ nn,j, (0,0), c'} ■ 

D.3 two-next-nearest-neighboring-triplet states (S = j, j = 0,1,2) 

When S = j' (j' = 0, 2) and p = 0, the trial states with two-next-nearest-neighboring triplets 

are 

I *nnnj',(o,o),a'> = I (0, 0)Ar/>+ | (0, 0)Ayj')+ | (0, 0)Bxj')+ | (0, 0)By/) (I+even), 

I ^nnn,j>,(o,o),b>) = \(0,0)Axj')+\(0,0)Ayj')-\(0,0)Bxj')-\(0,0)Byf) (I + even), 

I *nnnj',(o,o),c) = I (0, 0)Ar/>- | (0, Q)Ayj')+ I (0, Q)Bxj')- I (0, 0)By/) (J+odd), 
I *nnnj',(o ) o),d'> = I (0, Q)Axj')- | (0, Q)Ayj')- I (0, Q)Bxj')+ I (0, 0)£yj') 

When 5=1 and p = 7r they are 

I ^nnn,i,(n,7r),a") = I (tt, tt) Arl)+ | (iT,ir)Ayl) + i{\ (tt,it)Bx1)+ I (7r,vr)5yl)} (I + even), 

\^nnn,i,(^),b") = I (tt.ttJAcI)- | (tt, 7r)Ayl) - i{| (7r,7r)Sxl)+ I (7r,7r)Syl)} (/+odd), 

I *nnn,i,(7r,7r),c") = I (tt, /t)Ax1)+ | (7r,7r)Ayl) + i{| (tt, 7r)Sxl) - | (tt, 7r)5yl) } (I+even), 

I *nnn,i,( w ,ir),ci"> = I (ir,ir)Axl)- | (tt, 7r)Ayl) - i{| (tt, 7r)£xl) - | (7r,7r)£yl)} (I + odd). 

Otherwise they are 

I ^nnn. j,(p,p),a) = I (p,p)^j)+ I (p,p)Ayj) (I+even), 
I ®nnn,j,(p,p),b) = I | {p,p)Ayj) (I+odd), 

j,(p,p),c) = I {P-,P)Bxj)+ | (p,p)Byj) (I + even), 
I *nnnj,(p,p),<i} = I (P,P)Bxj)- \ (p,p)Byj) (I + odd). 

Note that E nnn ji^ n ^ a = E nnn j,^ n ^^ d , E nnn ji^ n ^ b = E nnn ji^ n ^ c , where j' = 0,2, and 

E nnn ,l,(0,0),a = ^nnn,l,(0,0),d) ^nnn,l,(0,0),6 = ^rmn,l,(0,0),c> ^nnn,l,(ir,7r),a" = ^nnn,l,(7r,7r),c" an d 
1 — 



7/ . 



nnn,l,(ir,Tr),b" • C/ nran,l,(7r,7r),<i- 

The following should be kept in mind, too ( j' = 0,2). 

^- I *nnn,j',(0,0),o') = + I *nnnj',(0,0),a') ) I *nnnj',(0,0),a') = + I *nnnj',(0,0),a') ■ 
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I. 


-1 


^nnn,j>, (0,0), 6') z 


= + 1 


^nnn,j',(0,0),b')i 


v x 


^nnn,j>,(0,0),b>) = _ 


^nnn,j', (0,0), b')> 


I. 


1 


^rmn,j',(0,0),c') = 




^ nnn,f, (0,0), c')i 


v x 


^nnn,j',(0,0),c') = + 


*nnnj',(0,0),c')i 


I. 


- 1 


^nnn,j>, (0,0), d'} = 


= " 1 


^nnn,j',(0,0),d')i 


v x 


^nnn,j>, (0,0), d') = — 


*nn7ij',(0,0),d')) 




J_ 


*rmn,l,(0,0),a) 




*rmn,l,(0,0),a)> 


v x 


*nnn,l,(0,0),a) = — 


*rarm,l,(0,0),d)j 




L 


- *nnn,l,(0,0),b) 


= + 


*nnn,l,(0,0),b); 


v x 


*nnn,l,(0,0),b) = ~~ 


^nnTiil, (0,0), c) i 




L 


- 1 ^ ' nnn, .1,(0,0), c) 




1 ^ nnn, 1,(0,0), cl i 


v x 


^nnTijl^OjO)^) — 


*rmn,l,(0,0),&)j 




/_ 


1 *jlTlTl,l,(0,0),d) 


= + 


1 ^ nnn, 1,(0,0), d) i 


v x 


^nnn, l,(0,0),d) = _ 


^rann, 1,(0,0), a) • 



Appendix E: Trial states in the p y = case 

We limit ourselves to S = 1 and p x / here. The trial states are as follows. 

E.l one-triplet states 

I *one,l,(p,0),o> = I (P,0)A) +exp(^) | (p,0)B), 
I * O ne,l,(p,0),b) = I (P, 0)A) - exp(|) I (p, 0)5). 

Note that 

V s {| (p,0)A> ±exp(|) | (p,0)B>} = ±exp(-|){| (p, 0)A) ± exp(|) | (p,0)B)} . 
E.2 two-nearest-neighboring-triplet states 



^ ' nn,l,(p,0),a) = 


| (p,0)AR) + eM l {) 


1 (P,0)BR), 


^ nn,l,(p,0),b) = 


| ( p ,0)Ai?)-exp(|) 


1 (p, 0)5/2), 


^nrijl^PjO)^) = 


| (p,0)AL>+exp(|) 


(P,0)BL), 


^nnjl^PjO),^) = 


| ( p ,0)AL)-exp(|) 


(p,0)5L). 



Note that 

V y {\ (p,0)AZ) ±exp(|) | (p,0)BZ)} = ±exp(-|){| (p, 0)AZ) ± exp(|) | (p,0)BZ)}, 

with Z = L or Z = R. 

Also note that -E n n,l,(p,0),a = -Enn,l,(p,0),c and £ nn ,l,(p,0),b = ^nn,l,(p,0),<2- 
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trial state | S i+ J_ U (E m -E )/J 



*ran,0,(0,0),a') 











0.49 


*nn,0,(0,0),6') 









+ 


0.40 


*nn,0,(0,0),c') 







+ 


* 


0.93 


*nn,0,(0,0),d') 





+ 




* 


0.93 


*nnn,0,(0,0),a') 





+ 


+ 


+ 


0.80 


*nnn,0,(0,0),b') 





+ 


+ 




0.91 



*nn,l,(0,0),a') 

*nn,l,(°> )> 6 ') 
*nn,l,(0,0),c') 
*nn,l,(0,0),d') 



1 - - 

1 - - 

1 - + 

1 + - 



0.89 

+ 1.04 

* 0.40 

* 0.40 



Table I. The lowest energy for the trial states with different quantum numbers at the coupling ratio J'/J = 0.65. 
Asterisks (*) in the table denote that the state is not an eigenstate of the operator U. 
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Fig. 1. A schematic view of the Shastry-Sutherland model. Each site is occupied by a spin. The intra-dimer coupling 
is denoted by dotted lines (A-sublattice) or dot-dashed lines (B-sublattice) and the inter-dimer coupling by solid 
lines. 
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-0.5 



Fig. 2. Difference between Eo = — |iVJ (the energy of the singlet-dimer state) and the lowest energy obtained from 
the trial states with the momentum (0,0). The dotted line is the spin gap estimated by the perturbation theory 
up to the fifth order, 13) A/J = 1 - (J'/J) 2 - |(J'/J) 3 - §(J' '/ 'J) 4 + 4(J'/J) 5 . 
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Fig. 3. Dispersion of the energy difference when the coupling ratio J' /J — 0.5. 
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Fig. 4. Dispersion of the energy difference when J'/J = 0.6. 
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Fig. 5. Dispersion of the energy difference when J'/J = 0.65. 



